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ABSTRACT 

We study the lensing convergence power spectrum and its covariance for a standard ACDM cos- 
mology. We run 400 cosmological A^-body simulations and use the outputs to perform a total of 
1000 independent ray-tracing simulations. We c ompare the si mulat ion results with analytic model 
predictions. The semi-analytic model based on ISmith et al.l ([20031) fitting formula underestimates 
the convergence power by ^ 30% at arc-minute angular scales. For the convergence power spectrum 
covariance, the halo model reproduces the simulation results remarkably well over a wide range of 
angular scales and source redshifts. The dominant contribution at small angular scales comes from 
the sample variance due to the number fluctuations of halos in a finite survey volume. The signal-to- 
noise ratio for the convergence power spectrum is degraded by the non-Gaussian covariances by up to 
a factor 5 for a weak lensing survey to ~ 1. The probability distribution of the convergence power 
spectrum estimators, among the realizations, is well approximated by a ^^-distribution with broad- 
ened variance given by the non-Gaussian covariance, but has a larger positive tail. The skewness and 
kurtosis have non-negligible values especially for a shallow survey. We argue that a prior knowledge 
on the full distribution may be needed to obtain an unbiased estimate on the ensemble averaged band 
power at each angular scale from a finite volume survey. 

Subject headings: gravitational lensing - large-scale structure of the Universe - cosmology: theory - 
methods: numerical 
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1. INTRODUCTION 

Weak gravitational lensing provides a unique probe 
of the mass distribution in the Universe. By detecting 
the so-called "cosmic shear", one can directly measure 
the power spectrum of mass fluctuations on cosmolog- 
ical scales of tens or even hundreds of mega-parsecs. 
The significant detections of cosmic shear signals has 
been reported by variou s groups, (e.g. iBacon et al] 
2000t iKaiser et al.l [2OOOI: iVan Waerbeke et al.l 120001: 



Wittman et all 1200 0: Hama na et al.l l2003l an"dits abil- 



ity t o constrain cosmological parameters has been shown 
(e^ lJarvis et all l2006HSemboloni et al.l 120061 : iFu et all 

mm . 

Weak lensing can also be a powerful probe of the na- 
ture of dark energy. The growth rate of mass fluct ua- 
tions can be rneasured by "lensing tomogr aphy" (e.g.lHul 
IT999t IHutereill200l iTakada fc ,Ta"g[200l which in turn 
provides tight constraints on the dark energy equation 
of state. A number of wide-field surveys are planned 
for this purpose, suc h as Subaru Weak Lens Survey 
(|Mivazaki et al.ir2006( ). the Panoramic Survey Telescope 
& Rapid Response System (Pan-STARRS^), the Dark 
Energy Survey (DES^), the Large Synoptic Survey Tele- 
scope (LSST^), and the Joint Dark Energy Mission 
(JDEM^). 



Electronic address: masanori@a.phys.nagoy a-u.ac.jpl 
^ http:/ /pan-Starrs. ifa.hawaii.edu/public/ 
^ http:/ /www. darkenergysurvey.org/ 
^ http://www.lsst.org/ 
* http://jdcm.gsfc.nasa.gov/ 



However, to attain the full potential of future surveys, 
it is of great importance to employ adequate statisti- 
cal measures of weak lensing for estimating cosmological 
parameters, properly taking into account correlations of 
the observables between different angular scales and red- 
shifts, i.e. the covariances. Since most useful cosmologi- 
cal information in weak lensing is inherent in small angu- 
lar scales that are affected by nonlinear cluste ring regime, 
the co variance is by nature non-Gaussian. Ilchiki et al.l 
()2009| ) argued that the use of an inaccurate covariance 
matrix may result in a systematically biased parameter 
estimate. However, a modeling of the covariance matrix 
requires an accurate modeling of nonlinear structure for- 
mation, which is not so easy. 

There are several approaches to modelling the co- 
variance matrix for a given set of cosmological and 
survey parameters. The most accurate way to ob- 
tain predictions for weak-lensing surveys is to per- 
form ray-tracing simulations through large- volume, high- 
resolution iV-body s i mulat i ons of cosmic structure 
formation (iJain et al.l 120001: iHamana fc Meiliei) 1200X1: 



Vale fc Whitd 120031: IWh ite fc Vald 120041: iFosalba et al.l 
2008; Tcyssicr eFaLllloOg; Hilbcrt et al. 2009). Over a 
range of angular scales of interest, the non-Gaussian ef- 
fects can be signifi c ant in weak lensing measurements 
(IWhite fc Hd [20001 : ISemboloni et~al[ [2007t lEifler et all 
I2008D . However, rather expensive calculations are needed 
to create many independent realizations in this method. 
A less expensive way would be to use semi-analytic mod- 
els that are based on, for instance, the so-called "halo 
model approach" ([Coorav fc Hul 120011 : iTakada fc Jain! 
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[2009I) . 

The purpose of this paper is to study the convergence 
power spectrum and its covariance using ray-tracing sim- 
ulations for a concordance ACDM cosmology. In order to 
obtain an accurate covariance matrix, we perform 1000 
independent ray-tracing simulations. We compare our 
simulation results with the halo model predictions for 
the covariance matrix. We also study the cumulative 
signal-to-noise ratio for measuring the convergence power 
spectrum taking into account the non-Gaussian errors. 
Recently, Lcc & Pen (2008) studied the angular power 
spectrum of the SDSS galaxy distribution and showed 
that the signal-to-noise ratio integrated over a range of 
multipoles is two orders of magnitude smaller than the 
case of Gaussian fluctuations. Motivated by their find- 
ing, we examine how the cumulative signal-to-noise ratio 
for the lensing power spectrum measurement is degraded 
by non-Gaussian errors using our large number of simu- 
lation realizations. Furthermore, we will study how the 
convergence power spectrum estimates are distributed in 
different realizations: we will study the probability distri- 
bution of convergence power spectrum and then compute 
the higher-order moments, skewness and kurtosis. 

The structure of this paper is as follows. In Section [5] 
we describe the details of iV-body simulations and ray- 
tracing simulations. In Section [31 after defining the lens- 
ing power spectrum, we show the simulation results for 
power spectrum estimation and then compare the results 
with the analytical prediction. In Section [3] we study the 
power spectrum covariance using the simulations and the 
halo model, and also estimate the expected signal-to- 
noise ratio for the power spectrum measurement taking 
into account the non-Gaussian errors. In Section [5] we 
study the probability distribution of power spectrum es- 
timators in our simulations. Section [S] is devoted to con- 
clusion and discussion. 

2. NUMERICAL METHODS 
2.1. The Cosmological N-body Simulations 

We u se the parallel Tree-Particle Mesh code Gadget-2 
(jSpring el 2005) in its fuU Tree-PM mode. We employ 
256^ particles for each of two different simulation vol- 
umes, 240 and 480/i~^Mpc on a side, which are used for 
ray-tracing simulations for redshift ranges of z = [0, 1] 
and [1,3], respectively (see Fig. [Ij. The smaller vol- 
ume simulation for lower redshifts is used in order to 
have higher mass and spatial resolutions because nonlin- 
ear clustering is more evolving at lower redshifts. We 
generate the initial conditions following the standard 
Zel'dovich approximation. In this step we employed the 
linear matter tran sfer function computed from CAMB 
(jLewis et al.ll200 (?) . The initial redshift is set to Zinit — 40 
and 50 for the large- and small-box simulations, respec- 
tively. We perform 200 realizations for each sets and 
hence have a total of 400 realizations. 

We adopt the concordance ACDM model with matter 
density = 0.238, baryon density ilb = 0.042, dark 
energy density f^A — 0.762 with equation of state pa- 
rameter w = —1, spectral index rig = 0.958, the variance 
of the density fluctuation in a sphere of radius 8 h~^Mpc 
(78 = 0.76, and Hubble parameter h = 0.732. These pa- 
rameters are consistent with the WMAP 3-year results 
(|Spergel et al..,2007,l . 



We have checked that our simulation result for the con- 
vergence power spectrum agrees with the result using 
512 particles within 5 % at / ^ 10''. This is sufficient 
for our purpose, which is to study the power spectrum 
and the covariance down to arcminute scales. 

2.2. Ray-tracing Simulations 

Fig. [1] shows the design of our ray-tracing simulations. 
We place the small- and large-volume simulations to 
cover a light cone of angular size 5° x 5°, from redshift 
z = to z ~ 3.5, using the tiling technique developed in 
IWhite fc Hi] dlOOO) and Hamana & Mellier (2001). 

We use the standard multiple lens plane algorithm 
in order to simulate the distortion and magnification 
of background light rays by foreground matter. Let us 
briefly describ e the procedure to trace light rays through 
A^-body data (jJain et al. 2000). In the standard multiple 
lens plane algorithm, the distance between observer and 
source is divided into A^ intervals, separated by comov- 
ing distance Ax- We adopt a flxed interval between lens 
planes by Ax = 120/i~'^Mpc (for this choice, the simula- 
tion box side lengths become multiples of 120/i~'^Mpc). 
The projected density contrast of the p-th plane is given 

by 

/•Xp-i+Ax 

Sp(0)= / dx6{xe,x), (1) 

where 5 is the three-dimensional density fluctuation flcld 
along the line-of-sight, 5 = p/ p — 1, and Xp-i denotes 
the x-position of {p-1 )-th lens plane, and is the two- 
dimensional vector denoting the angular position on the 
sky. The projected density fleld is computed on 4096^ 
grids by projecting A^-body particle distribution onto the 
lens plane based on the triangul ar-shaped cloud (TSC) 
assignment scheme (jHocknev fc Eastwood 1 988, ). An ef- 
fective two-dimensional gravitational potential of the p- 
th plane is related to projected density contrast via the 
two-dimensional Poisson equation 

VVp = 3(^^^ anSp. (2) 

This equation is solved to compute using the fast 
Fourier transform method making use of the periodic 
boundary conditions. Then, the flrst and second deriva- 
tives of ifjp are evaluated on each grid point. Next, 2048^ 
rays are traced backward from the observer point. The 
initial ray directions are set on 2048^ grids, which corre- 
spond to angular grid size of 5°/2048 ~ 0.15 arcmin. For 
each ray, we flrst computed ray positions on all the lens 
planes using the lens equation: 

^^^_gMXn^ (3) 
^ fK iXn) 

where fxix) is the comoving angular diameter distance 
given as function of x; /if(x) = X for ^ ^cit universe. 
The first and second derivatives of ipp on a ray position 
are linearly interpolated from four nearest grids on which 
they were pre-computed. The evolution equation of the 
Jacobian matrix, which describes deformation of an in- 
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Fig. 1. — The design of our ray- tracing simulations. Dashed lines 
show the field-of-view spanning ±2.5 degree. 



finitesimal light ray bundle, is written as 
//f(Xp)/x(Xn - Xp) 



n-l 

p=l 



/k(x«) 



UpAp, (4) 



where I is the identity matrix, and Up is the shear tensor 
on the p-th lens plane defined by 



dxidxi 



(5) 



where Xi = x^i s-nd so on. The Jacobian matrix is usually 
decomposed as 



A = 



1 - K - 7i 
-72 + w 



-72 - w 

- K + 71 



(6) 



7!)"'^/^ is the magni- 



where k is convergence, I7I = (71 
tude of the shear, and w is a net beam rotation. Finally, 
the summation in Eq. (j4]) yields the Jacobian matrix. 
The light-ray positions and four components of the Jaco- 
bian matrix on desired source planes are stored. Source 
redshifts we consider in this paper are summarized in 
TableU] 

We perform 1000 ray-tracing realizations of the under- 
lying density field by randomly shifting the simulation 
boxes assuming periodic boundary conditions. In doing 
this, each simulation output is shifted in the same way 
to make several lens planes to maintain the clustering 
pattern of mass distribution within the simulation box 
Importantly, however, in order to have independent real- 
izations, wc do not use the same simulation output when 
making each ray-tracing realization. Note that we use 
only one projection axis to all ray-tracing realizations, 
unlike many previous works in which three orthogonal 
directions are used to increase the number of realizations. 



^ Note that, for the separation between lens planes Ax 
120/i~^Mpc, we make two (four) lens planes from one simulation 
output for the redshift range z = [0, 1] (z = [1, 3] ) be cause we use 
simulations of different volumes as described in § 12.11 



TABLE 1 

Source redshifts employed for our ray-tracing 
simulations. 





exact values 


0.6 


0.588542 


0.8 


0.810822 


1.0 


0.996884 


1.5 


1.51902 


2.0 


1.99765 


3.0 


3.05725 



3. POWER SPECTRUM 
The mass density power spectrum Ps{k) is defined as 

{5{k)5*{k')) = i27TfSl{k - k')Ps{k). (7) 

Likewise to Eq. ([7]), one can define the convergence power 
spectrum Pk{1) as 

{k{l)n*{n)^{2nrSl{l-l')P^{l). (8) 

The conversion from the 3-D wave vector k to the 2-D 
angular wave vector I is done by the line of sight inte- 
gration using the Limber approximation (|LimbeH 119541 : 
IKaiseri [l99l . By using the Limber approximation, 
the convergence power spectrum is given by (see, e.g. 
'Barte lmann fc Schneidei1l2001t IVan Waerbeke fc Mellieil 
,2003.1 



Mx) 



;x 



(9) 



where xh is the horizon distance, defined as the comov- 
ing distance obtained for infinite redshift. The weight 
function W^(x) is now 



2c2a(x)- 



XH 



Ik{x) / dx'G(x') 



fK{x'-x) 
fKix') 



(10) 

where G(x) is the probability distribution of sources in 
comoving distance. For simplicity, all sources are as- 
sumed to be located at the same redshift z,, so that 



G{z) = 5d{,z - Zs). 



(11) 



The binned convergence power spectrum can be esti- 
mated from each realization as 



1 



l:\l\el 



(12) 



where the summation runs over modes whose lengths lie 
in the range I — Al/2 < < I + Al/2 for the assumed 
bin width Al. Throughout this paper we employ the bin 
width Alnl = 0.3. The number of modes around a bin I 
is approximately given by 



Ni^A,- a/(27r)^ 



(13) 



where As is the area of the two-dimensional shell around 
the bin I and can be given as ~ 2ttIAI ~\- 7r(AZ)2, and 
Vis is the survey area. Taking the average of Eq. over 
a number of realizations is expected to give the ensemble 

average expectation Pk{1) — {Pk{1))- 

Fig. [2] shows the convergence power spectrum obtained 
from our ray-tracing simulations for Zs = 1.0. We 
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Fig. 2. — Convergence power spectrum for sources at redshift 
Zs = 1.0. The result from our 1000 ray-tracing simulations is shown 
as the cross symbols with error bars (ztlcr variance). We also show 
the semi-analytic prediction computed from Eq. ^ using HaloFit 
to model the 3D mass power spectrum. 




k(h/Mpc) 

Fig. 3. — Top panel : Comparison of the dimensionless 3D mass 
power spectrum, A'^(fc) = P{k) /2tt'^ , obtained from 200 realiza- 
tions to the HaloFit prediction at 2 = and 0.92. Bottom panel 
: The fractional difference is shown in percent. The arrow shows 
the Nyquist wavenumber. 



compare i t with the semi- an alytic prediction computed 
using the ISmith et al.l (|2003D fitting formula (hereafter 
HaloFit) to compute the nonlinear matter power spec- 
trum Ps{k). The HaloFit underpredicts the convergence 
power at intermediate and small sca les, I > 300 0. A simi- 
lar disagreement is also found in Hil bert et al.! ( 2009) us- 
ing Millennium Simulation (,Springel et al.l 12005 ) which 
has higher mass and spatial resolutions than ours. 

To further explore the cause of this discrepancy, we 
study the three-dimensional mass power spectra using 
our iV-body simulations. Fig. [3] compares the power 
spectra obtained from 200 realizations with the HaloFit 
predictions at z = and 0.92. The arrow shows the 
Nyquist wavenumber. The HaloFit results are approxi- 
mately 5 — 10% l ower t han the results from simulations. 
iHeitmann et al.l (|2008f ) also report a similar disagree- 
ment. It appears that the discrepancy in the convergence 
power at high I is owing to inaccuracy in HaloFit. 

We also examine the effect of smoothing used in ray- 
tracing simulations by varying the grid size. Fig. |4] com- 
pares the convergence power spectra calculated for two 
different grid sizes, 2048^ and 4096^ grids, respectively. 
Clearly, the coarser grid size yields a smaller power at 
/ > 3000. In the bottom panel, we show that the differ- 
ence reaches ~ 5% at I ^ 3000. When we compute the 
projected density field, we would naively expect that a 
finer grid size provides a higher resolution in the lens- 
ing convergence map, if the original A^-body simulation 
has a sufficient spatial resolution. Further halving the 
grid size, i.e. changing the grid number to 8192^ from 
4096^, is similarly thought to give a better accuracy up 
to higher multipoles. We have checked the difference 
reaches ^ 5% at I ~ 6000. Therefore, we conclude 
that estimating power spectrum from the projected den- 
sity fields on 4096^ grids is sufficiently accurate up to 
I ~ 6000. In the following section, we focus on the power 
spectrum information up to ^ '--^ 6000. 

4. COVARIANCE MATRIX 
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10= 103 IQ4 IQ5 

multipole 1 

Fig. 4. — The effect of grid size in ray-tracing simulation on the 
power spectrum estimation. Top panel : The two spectra computed 
from 1000 realizations employing 2048^ and 4096^ grids. Bottom 
panel : The fractional difference is shown in percent. 



The covariance matrix of the convergence power spec- 
trum between Pk,{1) and Pk{1') is formally expressed as 
a su m of the Gaussian and non-Gaussian contri butions 
(e.g. iScoccimarro et al.lll999l : iCoorav &: Hiill2001h : 

cov[p«(o,Pk(0]^((^^«(o-^«(o) {p.in-p.in)) 

(14) 
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Fig. 5. — The halo model predictions for the ratio of the diag- 
onal non-Gaussian covariance components to the Gaussian errors 
as a function of multipole, where the Gaussian errors are given by 
Cov^ = P^COVC/sky'^Aln/). Note that we assume source redshift 
Zs = 1, the survey area f2s = 25 sq. degrees and the multipole bin 
width Aln/ = 0.3. In our halo model the non-Gaussian errors are 
given by the sum of t he d iffer ent halo terms of trispectrum con- 
tribution (see Egs. 1141 and I17I I and the halo sample variance (see 
Eg. I18I I. The thick solid curve shows the total contribution to the 
non-Gaussian errors, while the thin solid curves show each different 
term contributions as indicated by each label. 

where S^^, is the Kronecker deha function and T is the 
lensing trispectrum defined as 

(15) 

where we have introduced notation ^1234 — ii+i2+^3+'4- 
In the Limber approximation, T is a simple projection of 
the three-dimensional mass trispectrum Tg given by 



Jo 



dx 



X Ts 



I9 



I3 h 



fxix) fKix) fK{x) fKix) 



X 



(16) 



In Eq. (fTi|) . the first term describes the Gaussian er- 
ror contribution that has vanishing correlations between 
diS'erent multipole bins, whereas the second term de- 
scribes the non-Gaussian contribution arising from mode 
coupling due to nonlinear clustering. Both the terms 
scale with the survey area as oc l/fig- It should be also 
noted that the Gaussian term depends on the bin width 
AZ, whereas the non-Gaussian term does not (because 
J^(Pl/As « 1). Thus decreasing AZ increases the Gaus- 
sian contribution relative to the non-Gaussian errors. 

4.1. Halo Model Approach for the Covariance 

To make an analytic estimation of the lensing power 
covariance using Eq. (fTi)) . we need to model the mass 
trispectrum that can account for the non-linear cluster- 
ing at small angular scal es. In this paper we employ 
the halo model approa ch ()Seliakll200a iMa fc Frvl 120001 : 
IPeaco ck fc SmithI I200dl . also see Cooray & Sheth 2002 
for a thorough review). In the halo model, the power 
spectrum is given by a sum of two terms: the so-called 
1-halo term which describes correlation of dark matter 
particles within each halo, and the 2-halo term which 
describes correlation of particles in different halos. Like- 
wise, the trispectrum consists of four terms, from one to 



four halo terms (jCoorav fc Hull2001[ ): 



Ih 



rp2h 



3h 



rpih 

^5 ; 



(17) 

where we have suppressed the arguments (fei, ^2, ^a, '^4) 
of Ts ■ These four terms contribute to the power at each I 
differently. The 1-halo term gives dominant contribution 
in the nonlinear clustering regime. 

To complete the halo model approach, we need suitable 
models for the three ingredients: the halo mass density 
profile, the mass function of halos, and the biasing of halo 
distribution, each of which is specified by halo mass m 
and redshift z for a given c osmological model. For these 
we employ an NF W profile ( Navarro et al.llT997l ) , and the 
fitting formulae for mass function and halo bias derived 
in She th fc TormenI (|1999D . The details of our halo mode l 
implementati on can be found in PTakada fc JainI (|2003aD 
(also see Taka da fc j"ainll2003bD . 

However it turns out there is additional contribu- 
tion to the non-Gaussian covariance, which becomes 
significant on small scales as descr i bed b elow. As 
first pointed out in Hu fc Kravtsovl (|2003D . the halo 
clustering causes additional sampling variance to the 
covariance due to the statistical fluctuations in the 
number of halos sampled by a finite survey volume 
(also see [T akada fc Bridle 2007; Nevrinck et al. 2006|; 
IZhang fc Sh eth 2007; Reid & Spergel 2009). In simpler 
words, if massive halos happen to be more (less) in a 
surveyed region, the measured power spectra would very 
likely have greater (smaller) amplitudes than expected 
from the ensemble average. 

Accordin g to the formulation developed in Appendix in 
iTakada fc Bridle (2007) (also see the discussion around 
Eq. (7) in Zhang fc ShethI I2007t ). the sample variance 
to which we hereafter refer as the halo sample variance 
(HSV) is expressed as 

cov„,vR(o,p.(ni^f'<ix(|)n)' 



dn 



dM'^b{M')\kM'{l')? 



dM 



"kdk 



X —P^{k-x) W{kxQs) , (18) 



where d^V/dxdfl is the comoving volume per unit 
solid angle and unit radial comoving distance, given as 
d^V/dxdil — x^ for a flat universe, Xs is the comoving 
distance to a source redshift z^, dn/dM is the halo mass 
function, b{M) is the halo bias parameter, and «:m(0 is 
the angular Fourier-transform of the conv ergence field for 
a halo with mass M (see Eqs. 28 and 31 in lTakada fc JainI 
[200 3a, for the definition). The quantity Ps{k) is the lin- 
ear 3D mass power spectrum, and W(x) is the Fourier 
transform of the survey window function; for this we sim- 
ply employ W{kxQs) = Ji{kxQs) / {kx^s) {Ji[x) is the 
first-order Bessel function) assuming = xAWtt for a 
given survey area Og. Thus we have assumed that the 
survey area is sufficiently large and the number fluctu- 
ations of halos are in the linear regime. Note that this 
sampling variance contribution does not necessarily scale 
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with l//sky unlike other covariance terms. The sample 
variance depends on /sky via the shape of linear power 
spectrum. For a CDM spectrum it becomes smaller with 
increasing f2s, and decreases faster than the other covari- 
ance terms that have the scaling of fi"^ or f~^^ if fig is 
greater than a few hundreds square degrees, for multi- 
poles of interest (l > 1000). 

More exactly Eq. ^TE\i was derived by replacing the 
function S'(b) in Eq. (Bl) of Taka da fc Br idle (2007) with 
|km(';2)|^, which is the lensing power spectrum contri- 
bution due to a halo with mass M and at redshift z. In 
fact the first term on the r.h.s. of Eq. (Bl) corresponds 
to the 1-halo term of non-Gaussian errors in Eq. p^ . 
while the second term in Eq. (Bl) yields Eq. p^ . The 
contribution of Eq. (jlSp arises for any finite-volume sur- 
vey because the halo distribution has modulations due to 
the biased density fluctuations over the survey window. 
The full derivation of non-Gaussian covariance within the 
context of the halo model approach is beyond the scope 
of this paper and will be presented elsewhere. 

To obtain a more physical insight, it would be useful to 
note that the sample variance (fT8|) is roughly expressed 
as 

CovHsv - Paf^,ies)P^''{l)P^\n, (19) 

where b is the halo bias averaged over halo masses and 
redshift interval, crrms(0s) is the rms of angular mass den- 
sity fluctuations for the survey area, and P^^H) is the 
1-halo term of the convergence power spectrum. Here 
the combination of &(Trms(0s) gives the rms fluctuations 
in the number of halos in the survey area. Thus the 
sample variance strength is proportional to combined ef- 
fect of the convergence power spectrum and the number 
fluctuations of massive halos due to the large-scale mass 
density fluctuations. The sample variance (|18p is vanish- 
ing if the halo distribution does not have any clustering, 
i.e. is completely random (corresponding to the limit 
barms 0). Note that, on the other hand, the 1-halo 
term of the trispectrum accounts for the Poisson contri- 
bution to the sample variance arising from the discrete- 
ness nature of halo distribution. 

Fig. O shows the halo model predictions for the di- 
agonal non-Gaussian covariances relative to the Gaus- 
sian errors as a function of multipoles. Note that we 
consider source redshift Zs = 1 as a representative ex- 
ample, but the results are very similar for other source 
redshifts we consider in this paper. The different halo 
term contributions to the covariance are more impor- 
tant on small multipoles, while the 1-halo term becomes 
increasingly significant with increasing multipoles. Im- 
portantly the sample variance contribution due to the 
number fluctuations of massive halos, given by Eq. (fT8|) . 
becomes dominant over other non-Gaussian errors at 
high multipoles / ^ 1000, boosting the non-Gaussian er- 
ror strengths by an order of magnitude up to Z ~ 10** 
from the estimate without this new effect. For high mul- 
tipole limit, the ratio of the diagonal non-Gaussian co- 
variance to the Gaussian errors has an asymptotic seal- 
ingasCovNG/Cov« ^ b^a^^^P^^OP^^ /[2P^^il)yNi] ~ 
b'^a^^,ies)Ni/2 oc PA\nl. Hence, Cov^^/Cov^ cx P for 
a constant bin width A In Z as implied from the results 
around / ~ 10^ in Fig. [5l because 6CTrms(0s) is constant 
for fixed survey area and cosmological model. 
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Fig. 6. — The diagonal components of the convergence power 
spectrum covariance for Zs = 0.6,0.8,1.0,1.5,2.0 and 3.0, respec- 
tively. The results are divi ded by the expected Gaussian covariance 
(the first term in Eq. II14| |). Therefore, the deviations from unity 
arise from the non-Gaussian errors. The symbols are the simula- 
tion results, while the solid curves are the halo model predictions. 
The shot noise contribution for source redshift Zs = 1.0 assuming 
Hg = 30 arcmin"^ and ctj = 0.22 for the mean number density and 
the rms intrinsic ellipticities, respectively. 

This contribution has been ignored in previous studies 
and we will below study the impact on the lensing power 
spectrum covariance comparing the halo model predic- 
tions with the simulation results. 

4.2. Diagonal Components of the Covariance Matrix 

Fig. [6] shows the diagonal components of the conver- 
gence power spectrum covariance as a function of mul- 
tipole. The values are divided by the expected Gaus- 
sian covariances for the power spectrum, which is esti- 
mated by inserting the power spectrum measured from 
the simulations into the first term on the r.h.s. of 
Eq. (jl4p . Therefore, the deviations from unity arise from 
the non-Gaussian error contribution. The different sym- 
bols are the simulation results for different redshifts. It 
is clear that the non-Gaussian errors start to be signif- 
icant at multipoles from a few hundreds, and the non- 
Gaussian errors are greater for lower source redshifts due 
to stronger nonlinearities in the large-scale structure. 

For comparison, the solid curves show the halo model 
predictions including the sample variance contribution 
due to the number fluctuations of halos in the simula- 
tion volume as described around Eq. ^5)1 . Note that, 
to obtain the halo model prediction, the survey area is 
set to fts = 25 deg^ as assumed for the ray-tracing sim- 
ulations. Rather unexpectedly the empirical halo model 
fairly well reproduces the simulation results over a wide 
range of multipoles and for redshifts we have considered. 
It should be noted that the agreement cannot be found 
if the sample variance (|18p is not included: the sam- 
ple variance is dominant over other non-Gaussian covari- 
ance terms at multipoles I ^ 1000 (sec Fig. [5]). The l- 
dependence of the ratio at these high multipoles is ap- 
proximately given as Cov^'^/Cov^ oc for all source 
redshifts, shallower than the asymptotic behavior cx 
discussed below Eq. p9|) due to the residual contribu- 
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Fig. 7. — The correlation coefficients rlP^il), Pk{1')] as a function of I for a given I', where I' is chosen to I' = 131 (left panel), 588 
(middle) and /' = 2635 (right), respectively. The solid curves denote the halo model predictions. Although the simulation and halo model 
results are in fairly good agreement, the simulations results display slightly greater correlation strengths for high multipolos and at lower 
redshifts. 
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Fig. 8. — Comparison of the covariance correlation matrices predicted by the halo model (upper triangular parts of the matrices) and 
those obtained from our 1000 ray-tracing simulations (lower triangular parts) for Zs = 0.6, 1.0,2.0 and 3.0. 
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tions of other non-Gaussian covariance terms. 

In practice, the shot noise contamination due to intrin- 
sic galaxy elhpticities contributes to the diagonal term of 
the covariance. The dashed line denotes the contribution 
for source redshift Zg = 1.0 assuming fig — 30 arcmin^^ 
and (Te = 0.22 for the mean number density and the rms 
intrinsic cUipticities, respectively, which are typical num- 
bers for a ground-based weak lensing survey such as the 
planned Subaru weak lensing survey. It is found that the 
shot noise becomes significant compared to the cosmolog- 
ical non-Gaussian contributions at very high multipoles, 
thereby making the covariance be closer to Gaussian in 
the multipole range. 

4.3. Off-Diagonal Components of the Covariance 
Matrix 

The correlation coefficients of the convergence power 
spectrum covariances quantify the relative strengths of 
the off-diagonal components to the diagonal components. 
The correlation coefficient is defined as 



r[F„(0,P.(/')] = 



Cov[P«(0,P.(n 



VCov[P,(/), P,(/)]Cov[P«(/'), Pk(^')] 

(20) 

Thus the the correlation coefficient is defined so as to 
give unity when I = I'. For off-diagonal components 
r — > 1 implies strong correlation between the two spectra 
of different multipoles, while r means no correlation. 

Fig. [7] shows the correlation coefficients r{l,l') as a 
function of / for a given I' and at Zs = 0.6, 1.0, 2.0 and 
3.0. Note that the results depend on the bin width. The 
solid curves denote the halo model predictions. The halo 
model predictions fairly well reproduce the simulation 
results. A closer look implies a sizable disagreement for 
very high multipoles and at lower redshifts. 

Fig. [5] shows the correlation coefficient matrices at 
Zs = 0.6,1.0,2.0 and 3.0 in a two-dimensional multi- 
pole space of (1,1'). The upper triangular parts of the 
matrices are the halo model prediction, while the lower 
triangular parts are the simulation results from our 1000 
realizations. The correlations are generally stronger at 
higher / and at lower redshift, as expected. 

4.4. Signal-to- Noise Ratio 

The obtained covariance matrices can be used to esti- 
mate the expected signal-to-noise {S/N) ratio for mea- 
suring the lensing power spectrum. The cumulative 
signal-to-no ise ratio can be defined (e.g. iTakada fc J aim 
I200I [2n?)9l ) as 



= J2 p.(.i)cov-'[p^{i),p^in]p,in, (21) 



Ll'<l„ 



where Cov~^ is the inverse of the covariance matrix and 
the power spectrum information over 72 ^ l < Imax is 
included (/ = 72 is the fundamental mode of our ray- 
tracing simulations. If ~ 27r/5° ~ 72). The signal-to- 
noise ratio is independent of the bin width, as long as 
the convergence power spectrum does not rapidly vary 
within bin width. 

Fig. [S] shows the S/N for the convergence power spec- 
trum as a function of maximum multipole ^max for z^ = 
0.6, 0.8, 1.0, 1.5, 2.0 and 3.0. The dotted line shows S/N 




maximum multipole: 1^^^^^ 

Fig. 9. — The cumulative signal-to-noise ratio (S/N) for the con- 
vergence power spectrum is shown as a function of maximum mul- 
tipole (max at 2s = 0.6,0.8,1.0,1.5,2.0 and 3.0, where the power 
spectrum information over a range of multipoles 72 < I < I 

max is 

included. The solid curves show the halo model predictions, while 
the dotted line is the result for the Gaussian covariance case. 

when only Gaussian errors are included. Our simulation 
results suggest that the S/N begins to deviate signifi- 
cantly from that of the Gaussian case. It increases slowly 
with increasing Imax in the non-linear regime. The S/N 
for low redshift surveys does not increase significantly 
at multipoles from a few hundreds to 1000. This im- 
plies that there is little gain in the S/N even if including 
modes at the larger I, as has been found in the pre- 
vious works (iRimes fc H amilton 2005; Hamilton et al 
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Taka hashi et al.ll2009l : (Pore et al.ll2009( l. The simulation 
results show that the S/N is degraded by non-Gaussian 
covariances by up to factor 5 for source redshift Zs — 1.0. 
iLee fc PenI ()2008f ) measured the Fisher information con- 
tent for the angular power spectrum of SDSS galaxies, 
which is equivalent to the S/N defined in Eq. pTjl . They 
found a similarly significant saturation of the S/N . 

The simulation results can be compared with the halo 
model predictions, where the non-Gaussian errors are 
computed from the halo model, and the power spec- 
trum and the Gaussian term of covariance are taken from 
the simulation results in the S/N calculation. The halo 
model predictions are again in a good agreement with the 
simulation results, and in particular well capture com- 
plex dependences of the S/N on multipoles and source 
redshifts. 

For planned weak lensing surveys most important 
range of multipoles are around I ~ 1000 in order to avoid 
complications due to effects of strong nonlinear clustering 
and gas dynamics on mass power spectrum. Over such a 
range of multipoles the cosmological non-Gaussian errors 
are dominant over the shot noise due to intrinsic galaxy 
shapes, therefore the effect needs to be properly taken 
into account to obtain unbiased, secure extraction of cos- 
mol ogical parameters frorn the measured power spectrum 
(see lTakada fc JainI (|2009f l for the similar discussion). 

5. PROBABILITY DISTRIBUTION OF THE 
CONVERGENCE POWER SPECTRUM ESTIMATOR 
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Fig. 10. — Probability distribution of the convergence power spectrum estimators among the 1000 reaUzations for Zs = 0.6 (left panel) 
and 3.0 (right), respectively. The solid and dashed curves show the x^- and Gaussian-distributions with zero mean, respectively, where the 
variance is set to the to the diagonal covariance components measured from the simulations. 



In this section, we study the probability distribution of 
the convergence power spectrum estimator in order to 
see how the non-linear growth causes a non-Gaussian dis- 
tribution in the convergence power spectrum estimators. 
It should be noted that the convergence power spectrum 
covariance simply reflects the width of the full distribu- 
tion at each I. 

Fig. [TO] shows the probability distribution of the con- 
vergence power spectrum estimators among 1000 re- 
alizations. We measure the probability distribution for 
{Ni/2y/'^{Pf^/PK — 1) so that the mean and variance of 
the distribution are equals to zero and unity in the linear 
regime or if the convergence field is Gaussian. The prob- 
ability distribution is normalized so as to give unity when 
integrated over all P«; values. For large I values, the dis- 
tribution is broadened due to non-linear evolution. The 
solid and dashed curves show the x^-distribution and the 
expected Gaussian distribution, respectively, where the 
variance for each of the distributions is set to the di- 
agonal term of covariance measured from simulations at 
each /, i.e. the variance includes the non-Gaussian co- 
variance contribution as given in Fig. [6l To be more 
explicit, assuming that the estimate Pk obeys the x^- 
distribution, the mean and variance are set to Pk(0 
and Cov[PyO,^«(0] (replacing P(fc) ^ P^jl), Nk/2 
P^Cov[P«(0,Pk(0] in Eq. (Bl) in ITakahashi"inil 
(liooa i). Fig. [TO] shows that the probability distribu- 
tion is well approximated by the x^-distribution, but 
display a larger positive tail than expected from the x^- 
distribution. One can see that the estimators have a 
skewed distribution, even for the low multipoles such as 
I ^ 100 where the lensing fields are more in the linear 
regime. 

The deviations from a Gaussian distribution can be 
quantified by studying skewness 5*3 and kurtosis 5*4 of 



the distribution in Fig. [TO] 

((A(0-^«(Of) 



53 



^4 = 



p.iwr'^ 



{{p.{i)-p.{i)W 



- 3. 



(22) 



If the convergence field is a Gaussian random field, which 
is a good approximation in linear regime, the convergence 
power spectrum estimator P^ of a given I exactly obeys 
the x^-distribution. In this case, the skewness and kur- 
tosis can be analytically computed as 



^3 



/4Cov[PyO,^«(0] 



Sa 



6 Cov[Py 0,^.(0] 



(23) 

,-1/2 



Note that S3 and 5*4 scale with survey area as S3 oc Qs 
and Si tx as Gov oc ft~^. 

Fig. [TT] shows the simulation results for ^3 and S4 as a 
function of multipole I. Note that the results in Fig. [TT] 
are for a survey area of fig = 25 degree^. The solid curves 
are the theoretical predictions of Eq. (|23|) . The model 
well describes the simulation results for Zs = 3.0, while 
the results for Zs — 0.6 show significant non-Gaussian cu- 
mulants over a range of multipoles we have considered, 
due to stronger non-linearities. For Zs — 3.0, both S3 
and 5*4 asymptote to zero at high I, i.e. the probabil- 
ity distribution approaches to a Gaussian distribution at 
high / due to the central limit theorem. 

Since the skewness and the kurtosis have non-negligible 
values at multipoles relevant for future surveys, a prior 
knowledge on the full distribution may be needed to ob- 
tain an unbiased estimate on the ensemble averaged band 
power at each I bin. 

6. CONCLUSION AND DISCUSSION 

Accurate statistics are essential in the likelihood anal- 
ysis for future precision cosmology. We can exploit the 
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Fig. 11. — The skewness (top panel) and the kurtosis (bottom) 
of the convergence power spectrum distribution as a function of 
multipole I at 2s=0.6 and 3.0. The solid curves are the theoretical 
predictions from Eq. (I23I I. 

full potential of upcoming high quality data, only if we 
use appropriate statistical methods. For weak-lensing 
surveys, non-linear gravitational evolution of large-scale 
structure can significantly compromise cosmological pa- 
rameter estimations and thus needs to be modelled ac- 
curately. 

We have used ray-tracing simulations in order to study 
how the non-Gaussian covariance varies with scales and 
redshifts for the standard ACDM cosmology. We have 
performed a total of 1000 independent ray-tracing simu- 
lations using 400 cosmological iV-body simulations. The 
non-Gaussian errors become more significant on smaller 
scales and at lower redshifts. The cumulative signal-to- 
noise ratio (S/N) for measuring lensing power spectrum 
is degraded due to non-Gaussian covariance by up to a 
factor of 5 for a weak lensing survey to Zs ~ 1. We show 
that the simulation results are fairly well described by 
the halo model prediction including additional contribu- 
tion due to the statistical fluctuations in the number of 
halos in a finite survey volume. 



We also study the probability distribution of the con- 
vergence power spectrum estimator among 1000 realiza- 
tions. The probability distribution has a large skewness 
especially for shallow surveys, which is likely due to non- 
linear gravitational evolution. Therefore, a prior knowl- 
edge on the full distribution may be needed to obtain an 
unbiased estimate on the ensemble averaged band power 
at each I. Overall, the non-Gaussian errors likely cause 
best-fitting parameters to be systematically biased, if the 
model fitting is done improperly assuming the Gaussian 
covariances. Therefore it is clearly needed to develop an 
appropriate method for parameter estimations from the 
measured power spectrum taking into account the non- 
Gaussian errors. 

The most conventionally used statistical measure is the 
cosmic shear correlation function. An invaluable feature 
of the correlation function is that it does not require 
non-trivial corrections for survey geometry and masking 
effects. Therefore, it is useful to estimate the covariance 
matrix of real space correlation function and to derive 
fitting formula to calibrate the full covariances for an 
arbitrary survey area. These issues will be studied in a 
subsequent paper. 
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